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LONG  TERM  GOALS 

The  long  term  goal  of  this  research  is  to  improve  the  accuracy  of  the  computed  multi-scale  flow 
physics  in  geometrically  and/or  hydrodynamically  complex  oceanic  and  coastal  ocean  environments 
through  refinements  in  the  defined  physics,  domain  definition  and  computational  grid  resolution.  The 
focus  is  in  particular  on  the  improved  coupling  of  wind  generated  short  wave  and  circulation  models 
within  the  framework  of  unstructured  adaptive  grid  models  for  oceanic  and  coastal  waters. 

OBJECTIVES 

The  objective  of  this  project  is  the  dynamic  coupling  of  the  ADCIRC  circulation  model  and  the  short 
wind  wave  model  SWAN  using  unstructured  computational  meshes.  In  recent  experience  with 
coupling  the  continuous  Galerkin  based  version  of  ADCIRC  with  short  wind  wave  models,  it  has  been 
found  that  wave  transformation  zones  require  high  levels  of  resolution  in  order  to  correctly  capture  the 
wave  radiation  stress  forcing  to  the  circulation  model.  However  the  wave  transformation  zone  tends  to 
shift  depending  on  the  direction  and  period  of  the  waves.  Thus  h-p  adaptive  discontinuous  Galerkin 
(DG)  based  solutions  are  being  developed  to  optimize  the  application  of  high  resolution  zones  to 
correctly  capture  the  wave  coupling  within  these  zones  without  over-resolving  adjacent  areas. 
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APPROACH 


We  have  focused  on  three  areas  to  accomplish  the  defined  objectives.  The  approach  consists  of  the 
following: 

•  The  first  area  is  implementing  and  evaluating  h-p  adaptive  schemes  for  Discontinuous  Galerkin 
(DG)  solutions  to  the  shallow  water  equations  (SWE)  in  anticipation  of  coupling  with  a  mesh 
adaptive  version  of  the  SWAN  wave  model  which  is  under  development.  This  will  improve  the 
unstructured  grid  capabilities  for  the  ADCIRC  circulation  model  by  implementing  /7/1-adaptive 
methods  to  ensure  that  sufficient  grid  resolution  is  provided  for  all  relevant  flow  scales.  The 
present  implementation  of  ADCIRC  is  based  on  both  CG  and  DG  finite  element  methods. 

Currently  we  are  intensively  developing  the  DG  based  algorithms.  DG  algorithms  are  particularly 
well  suited  for  both  propagation  and  advection  dominated  problems  with  or  without  sharp  gradients 
in  the  forcing  function,  bathymetry,  and/or  flow.  DG  methods  inherently  preserve  mass  perfectly 
on  an  elemental  level,  which  make  them  ideal  for  coupling  flow  and  transport  models.  Example 
problems  include  flows  with  strong  eddies  such  as  those  issuing  from  and  through  inlets;  flows 
through  rapidly  varying  bathymetry  such  as  canyons  or  deep  scour  holes  in  inlets;  and  flows  with 
sharp  fronts  such  as  tidal  bores,  ebb  tide  -  waves  interactions  in  inlets  and  wetting  fronts.  DG 
methods  are  conceptually  similar  to  finite  volume  methods  although  DG  methods  are  readily 
implemented  with  higher-order  bases.  DG  methods  are  also  ideal  since  non-conforming  h  and  p 
refinement  and  adaptivity  can  be  implemented. 

•  The  second  area  is  evaluating  the  influence  of  grid  resolution  in  the  SWE  circulation  model  in 
relation  to  the  wind  wave  model  in  the  computed  wave  radiation  stress  set  up.  In  addition,  we  are 
studying  the  sequencing  of  interpolation  and  the  wave  radiation  stress  gradient  computation  and  the 
influence  of  both  the  wave  and  surge  model  grids  in  this  context. 

•  The  third  area  is  to  improve  the  wave  -  circulation  model  coupling  which  was  bottle  necked  by 
our  single  processor  interface  for  communication  between  the  two  models.  We  are  designing  all 
paradigms  for  very  large  scale  parallelism  on  distributed  memory  machines,  anticipating  thousands 
of  processors.  We  are  migrating  all  interpolator  functions  to  the  individual  processors  and 
designing  direct  global  to  local  and  local  to  local  processor  communication  between  the  wave 
model  and  circulation  model  processors.  We  are  investigating  the  possibility  of  ADCIRC  and 
SWAN  running  on  identical  grids  on  the  same  processors,  entirely  eliminating  the  inter-model 
processor  to  processor  communication. 

WORK  COMPLETED 

A  DG  model  system  has  been  developed,  implemented,  and  tested  for  solving  the  two-dimensional 
shallow  water  equations,  sediment  transport,  and  passive  scalar  transport  equations  on  unstructured 
triangular  grids.  The  code  features  both  standard  h-  (grid)  and  p-  (polynomial  order) 
refmement/adaptivity  options  (the  current  model  is  adaptive  only  in  p).  The  models  have  been 
rigorously  verified  on  a  series  of  idealized  test  cases  using  both  the  /^-refinement/ adaptive  strategies, 
and  a  recently  parallelized  version  of  the  code  is  aiding  in  the  validation  process  of  the  model  on  full- 
scale  applications. 
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With  regard  to  some  specific  algorithmic  details  of  the  model,  through  the  use  of  an  orthogonal, 
hierarchical  basis  and  optimal  triangular  quadrature  rules,  an  efficient,  matrix-free  DG  method  was 
developed  and  implemented  as  the  unified  solution  technique  for  both  the  shallow  water  and  transport 
models.  The  hierarchical  nature  of  the  basis  that  was  employed  gives  the  models  a  flexibility  that 
allows  for  easy  implementation  of  ^-refinement  and  also  /7-adaptivity.  Additionally,  flexible 
programming  of  the  models  allows  for  the  use  of  arbitrary  stage  and  order  strong-stability-preserving 
(SSP)  Runge-Kutta  time  discretizations  to  be  used  in  conjunction  with  the  DG  spatial  discretization. 
This  allowed  for  a  detailed  study  of  the  computational  efficiency  of  a  number  of  different  time 
stepping  schemes. 

In  the  critical  area  of  inter-model  communications  we  have  also  made  substantial  progress.  Due  to  the 
bottlenecking  that  occurred  in  our  preliminary  implementation,  with  the  circulation  and  wave  models 
running  on  two  sets  of  processors  and  inter-model  communications  and  interpolation  handled  by  one 
master  processor,  we  have  developed  improved  strategies.  We  have  completed  implementation  of 
localized  interpolation  based  on  global  input  files  together  with  on-the-fly  globalization  by  the 
ADCIRC  model.  This  allows  for  testing  of  difficult  wind/wave  coupling  scenarios.  However  in 
anticipation  of  more  complex  scenarios  we  are  developing  processor  to  processor  communications 
between  the  models.  We  are  also  investigating  the  possibility  of  the  circulation  and  wave  models 
running  on  identical  grids  on  the  same  processors,  entirely  eliminating  the  inter-model  processor  to 
processor  communication. 

RESULTS 

The  following  results  are  noted  regarding  this  work: 

•  The  h  and  p  convergence  properties  of  the  hp-DG  circulation  model  were  demonstrated  for  both 
linear  and  highly  nonlinear  problems  that  included  the  effects  of  source/sink  terms,  spatially 
varying  bathymetry,  and  time  varying  boundary  conditions.  It  was  found  that  standard  h- 
refinement  for  a  fixed  p  leads  to  optimal  convergence  rates  of  p+ 1  in  both  surface  elevation  and 
velocity  (see  Figure  1),  while  exponential  convergence  rates  are  observed  for  /^-refinement  for  a 
fixed  h  (see  Figure  2). 

•  /7-refinement  versus  /7-refinement  was  investigated  for  a  simple  linear  problem  and  a  highly 
nonlinear  problem  with  complex  two-dimensional  flow  structure.  In  both  cases,  it  was  found  that 
significant  cost  efficiency  benefits  could  be  obtained  by  using  /7-refinement  instead  of  standard  h- 
refinement.  In  the  linear  case,  a  quantitative  analysis  comparing  model  run  times  to  error  levels 
demonstrated  convincingly  the  superiority  of  /7-refinement  for  smooth  problems.  This  fact  is 
shown  in  Figure  3,  which  plots  model  run  times  versus  error  levels  for  a  range  of  grids  and 
polynomial  orders,  where  it  can  clearly  be  seen  that  using  /7-refinement  is  more  computationally 
efficient  in  this  case.  In  the  case  of  the  nonlinear  problem  -  the  case  of  an  idealized  inlet  (see 
Figure  4)  -  a  qualitative  comparison  of  the  solutions  for  two  grid  refinements  and  a  range  of  p 
approximations  demonstrated  the  ability  of  /7-refinement  to  resolve  the  formation  and  advection  of 
strong  eddies  along  the  inlet  channel  even  on  very  coarse  grids.  This  is  shown  in  Figure  5. 
Higher-order  p  approximations  on  the  coarse  grid  gave  very  comparable  results  to  the  lower-order 
p  approximations  on  a  finer  grid,  the  latter  simulations  taking  two  to  four  times  longer  than  the 
former. 
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•  The  two-level  /(-adaptive  strategy  that  was  implemented  for  the  shallow  water  model  revealed 
the  ease  with  which  stable  and  accurate  p- adaptive  schemes  can  be  implemented  into  the  models, 
especially  given  the  choice  of  basis  functions  that  were  implemented.  Investigation  of  these 
schemes  in  the  context  of  the  idealized  inlet  problem  also  revealed  the  substantial  efficiency 
advantages  that  are  possible  with  these  techniques,  which  ran  in  36%  to  40%  of  the  CPU  time  as 
the  same  problem  using  the  higher-order  approximation  for  all  the  elements  for  all  time,  while 
obtaining  qualitatively  similar  results  (see  Figure  6).  Additionally,  this  paradigm  stands  to 
significantly  simplify  the  design  of  large-domain  unstructured  grids  that  can  adequately  resolve  the 
space  and  time  varying  flow  features  of  a  given  problem.  Starting  with  a  computational  grid  that 
has  been  designed  to  adequately  resolve  the  geometry  of  a  given  domain,  adaptivity  will  then 
automatically  adjust  the  elements  locally,  usually  based  on  some  local  characteristic  of  the  solution 
such  as  the  gradient,  in  order  to  provide  the  necessary  spatial  resolution  to  accurately  capture  the 
flow  field  as  it  evolves  in  space  and  time.  An  example  of  a  grid  evolving  is  shown  in  Figure  7, 
which  shows  snapshots  of  the  grid  for  the  idealized  inlet  at  different  times  during  a  simulation.  It 
can  be  noted  how  the  higher-order  elements  track  the  key  features  of  the  flow. 

•  In  a  one-dimensional  setting,  the  efficiency  of  a  number  of  SSP  Runge-Kutta  time 
discretizations  was  investigated.  The  class  of  SSP  methods  that  was  investigated  is  defined  by  the 
property  that  the  number  of  stages  s  is  greater  than  the  order  k  of  the  method.  From  analysis, 
conditions  on  the  time  step  (CFL  conditions)  for  the  linear  stability  of  the  methods  defined  using 
the  5  >  k  SSP  schemes  were  obtained  that  are  less  restrictive  than  those  of  the  “standard”  so-called 
RKDG  methods  that  use  s  =  k  SSP  Runge-Kutta  schemes.  The  improvement  in  the  CFL  conditions 
for  linear  stability  of  the  methods  more  than  offsets  the  additional  work  introduced  by  the  increased 
number  of  stages.  Thus,  more  efficient  RKDG  methods  that  possess  the  same  favorable  accuracy 
and  stability  properties  of  the  “standard”  RKDG  methods  were  obtained.  Numerical  results 
verified  the  CFL  conditions  for  stability  obtained  from  analysis  and  demonstrated  the  efficiency 
advantages  of  these  new  RKDG  methods.  The  theoretical  efficiency  advantages  that  are  possible 
with  these  schemes  over  the  standard  schemes  are  summarized  in  Table  1.  Numerical  results 
showed  excellent  agreement  with  these  theoretical  results. 

•  The  performance  of  the  models  was  also  demonstrated  and  evaluated  for  a  range  of  different 
problems  with  hydraulic  jumps  and  tidal  bores  (discontinuities  or  shocks)  and  advection  dominance 
in  both  one-  and  two-dimensions.  It  was  found  that  the  models  could  accurately  capture  the  steep 
fronts  and  discontinuities  in  surface  elevation,  and  also  concentration  or  bathymetry,  without  the 
introduction  of  spurious  oscillations,  as  is  common  when  using  standard  continuous  Galerkin  (CG) 
methods  for  these  types  of  problems. 

•  The  shallow  water  model  and  passive  transport  model  were  applied  to  a  full-scale  application  at 
Shinnecock  Inlet,  New  York.  The  models  were  shown  to  produce  stable,  non-oscillatory  results, 
and,  in  particular,  the  circulation  patterns  observed  in  the  vicinity  of  the  inlet  were  shown  to  be 
consistent  with  observational  data.  The  shallow  water  model  accurately  resolved  the  formation  and 
advection  of  eddies,  which  were  a  key  factor  in  determining  the  transport  and  fate  of  the 
contaminant  using  the  passive  transport  model  (see  Figure  8).  This  problem  is  a  good  candidate  for 
the  /(-adaptive  strategy  and  will  be  tested  in  the  near  future. 
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Open  Ocean  Boundary 


=  Land  Boundary 


Figure  4:  The  computational  domain  and  grid  of  the  idealized  inlet  test  case. 
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a.)  h {  ,p  =  \ 


b .)  \  ,p  =  2 


C •)  h  ,P  =  3 


Figure  5:  Computed  velocity  contours  and  vectors  for  the  idealized  inlet  test  case 
using  a  coarse  grid,  hi  with  polynomial  approximations  of  P  ~ l’  z’  ana  J 
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P  =  1: 


P  =  3: 


Adaptive  Method,  p  =  1  /  p  =  3 : 


Figure  6:  Comparison  of  the  P  ^,P  ^  and  adaptive 

solutions  for  the  idealized  inlet. 


p  =  1/ =  3 
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H  p  =  1 

p  =  3 

Figure  7:  Dynamic  evolution  of  the  elements  using  the  p-adaptive  algorithm. 
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Figure  8:  Transport  of  a  passive  quantity  at  Shinnecock  Inlet  at  various  times. 
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Table  1:  CFL  conditions  and  relative  efficiencies  of  the  SSP(s,  k)  RKDG  methods 

(optimal  schemes  are  boxed) 


Stages,  5 

SSP(s,2) 

RE  (%) 

SSP(s,3) 

RE  (%) 

SSP(s,4) 

RE  (%) 

CFL 

CFL 

CFL 

2 

0.3333 

— 

— 

— 

— 

— 

3 

0.5882 

15.00 

0.2000 

— 

— 

— 

4 

0.7611 

12.42 

0.3061 

12.88 

0.1696 

— 

5 

0.8966 

7.06 

0.4060 

17.90 

0.2152 

21.19 

6 

1.0089 

0.89 

0.4842 

17.39 

0.2747 

25.19 

7 

1.1052 

-5.09 

0.5667 

17.65 

0.3214 

26.12 

8 

1.1895 

-12.08 

0.6444 

17.24 

0.3707 

26.80 
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